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We present a microscopic derivation, without electrodynamical assumptions, of B{x,y,H{t)), 
M{H{t)), and Jc[H{t)), in agreement with experiments on strongly pinned superconductors, for 
a range of values of the density and strength of the pinning sites. We numerically solve the over- 
damped equations of motion of these flux-gradient-driven vortices which can be temporarily trapped 
at pinning centers. The field is increased (decreased) by the addition (removal) of flux lines at the 
sample boundary, and complete hysteresis loops can be achieved by using flux lines with opposite 
orientation. The pinning force per unit volume we obtain for strongly-pinned vortices, JcB ~ npfp '^, 
interpolates between the following two extreme situations: very strongly-pinned independent vor- 
tices, where J^B ~ Upfp, and the 2D Larkin-Ovchinikov collective-pinning theory for weakly-pinned 



straight vortices, where JcB 
pinning sites. 



Upfp. Here, Up and fp are the density and maximum force of the 



PACS numbers: 74.60.Ec, 74.60.Ge, 74.60.Jg 
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I. INTRODUCTION 

Flux distributions in type-II superconductors are com- 
monly inferred from magnetization and critical current 
measurements ^ and interpreted in the context of the 
Bean model or its variations. The Bean model, which 
has been widely used for over three decades, postulates 
that the current density in a hard superconductor (i.e., 
with strong pinning) can only have three values: — Jc, 
0, and -l-Jc, where Jc is the critical current density, 
which is independent of the local magnetic flux density 
B(a;,y,t). The Bean model and its many variants make 
no specific claims with regard to the microscopic mech- 
anism controlling the trapping of vortices. Bean's pos- 
tulate, Jc ^constant, was modified several times by Kim 
et al. §: J, ^ \/B [T]; Jc ~ l/(6o + B) [3^3^]; 
Jc ^ l/ibo + B + b2B^ + baB^ + . . .) [3^]; where b, are 
constants. On the other hand, Fietz et al. Q suggested 
that Jc exp{—B/bo) ; while Yasukochi et al. [|j sug- 
gested Jc ~ These, and other proposals made 
during the 1960s, were followed by several other phe- 
nomcnological modifications of Jc{H) during the follow- 
ing two decades [0J|]. A microscopic description, without 
assuming any particular i3-dependence of Jc, of these 
flux distributions — in terms of interacting vortices and 
pinning sites — can be very valuable for a better under- 
standing of commonly measured bulk quantities. 

One of the most effective methods of investigating the 
microscopic behaviour of flux in a hard superconductor is 
with computer simulations (see, e.g., |@,||, and references 
therein). In this paper, we present molecular dynamics 
(MD) simulations of the evolution of rigid flux lines in 
a hard superconductor. We first introduce our model 
for vortex-vortex and vortex-pin interactions as well as 



the corresponding antivortex interactions. We then in- 
vestigate the fiux profile which results from a varying 
applied field; from such fiux profiles we obtain full hys- 
teresis loops indicating that our model has the essen- 
tial microscopic ingredients underlying the experimen- 
tally measured macroscopic quantities. We also inves- 
tigate the behaviour of Jc{H) for a controlled range of 
pinning parameters. 



II. SIMULATION 

Our simulation geometry is that of an infinite slab of 
superconductor in a magnetic field applied parallel to the 
slab surface. Thus, demagnetization effects are unimpor- 
tant. We also treat the vortices as perfectly stiff, so that 
we need to model only a two-dimensional (2D) slice of the 
3D slab. Our system is periodic in the plane perpendicu- 
lar to the applied field, and we measure distances in units 
of the penetration length A. Here, we present results for 
a system of size 36A x 36A. The simulation, described in 
further detail below, consists of slowly ramping an exter- 
nal magnetic field. Flux lines enter the edge of the sample 
and their positions are allowed to evolve according to a 
T = MD algorithm. The resulting vortex distributions 
at any external field can then be deduced as a function 
of distance into the sample. 



A. Sample Geometry and Time-Dependent Field 

The actual sample region is heavily pinned, and ex- 
tends from position a; = 6A to a; = 30A (Fig. 1). Outside 
the sample itself is a region with no pinning which ex- 
tends from a; = OA to a; = 6A and from x = 30A to 



1 



X = 36A (with 36A = OA according to our periodic bound- 
ary conditions). This sample geometry is shown in the 
upper panels of Fig. 1 . Here, the sample (pinned) region 
occupies the central 2/3 of the system, and the unpinned 
region the outer 1/3. 

We simulate the ramping of an external field by the 
slow addition of flux lines to the outside unpinned re- 
gion. Because there is no pinning in this region, the flux 
lines there will attain a fairly uniform density, and we 
may define the applied field H as $o times this density. 
Flux lines from the external region will move into the 
sample through points at the sample edge where the lo- 
cal energy — as determined by the local pinning and vor- 
tex interaction — is low. Thus, our simulation models the 
real situation where vortices nucleate at such low-energy 
regions at the surface. 

Further, in a real superconductor, vortices near the 
surface are not expelled by their interior neighbors be- 
cause of a field-induced Meissner current flowing at the 
surface. Again, our external "bath" of vortices simulates 
this behavior by providing a balancing inward force, pro- 
portional to the external field, on those vortices near the 
sample boundary. 



B. Equations of Motion 



The force per unit length between two vortices lo- 



cated at Yi and r,- is 



r 



87r2A3 



A 



(1) 



We model the vortex- vortex force interaction in its exact 
form by using the modified Bessel function Ki. This 
force decreases exponentially at distances larger than 
A, and we cut off the (by then negligible) force at dis- 
tances greater than 6A. Further, we have cut off the 
logarithmic divergence of the force for distances less 
than 0.1 A. These cutoffs were found to produce negli- 
gible effects on the dynamics for the range of param- 
eters investigated. Thus, the force (per unit length) 
on vortex i due to other vortices (ignoring cutoffs) is 
= Y.'jZi fv Ki{\y, - Yj\l\) Vij . Here, the Yj are 
the positions of the N^, vortices within a radius 6A, 



(r^ - i"j)/|rj - r^l, = ±/o, and 



/o = 



87r2A3 



(2) 



The sign of the interaction is determined by /i,; we take 
— -K/o for repulsive vortex-vortex interactions and 
fy = — /o for attractive vortex-antivortex interactions. A 
vortex and antivortex annihilate and are removed from 
the system if they come within 0.3 A of one another [Q. 
Forces are measured in units of /o, lengths in units of A, 
and fields in units of ^o/A'^. 



We model the pinning potential Q as Np short-range 

(p) 

parabolic wells at positions . The equation of motion 
for a vortex moving with velocity v is f = rjv, where rj 
is the viscosity (w $o^c2//0ji, with p„ being the normal- 
state resistivity). Thus, the overall equation for the over- 
damped motion of a vortex subject to vortex-vortex and 
pinning forces is 



(3) 



where 



N. 



fc=i \ 



A 



(4) 



Here, 8 is the Heaviside step function, is the range 
of the pinning potential, and fp is the strength (maxi- 
mum pinning force) of each well, measured in units of 
/o. For all the simulations presented here £,p — 0.12A 
and 7] = I. The parameters we vary here are the pin- 
ning strength fp and the average distance between pin- 
ning sites dp (which determines the pinning density np 
via Up = 1/dp). Many other parameters can be var- 
ied, making the systematic study of this problem very 
complex. A more thorough investigation with different 
pinning-potential ranges, pinning potential-shapes, non- 
uniform strength distributions, and non-random pinning 
positions will be presented elsewhere. Here, the pinning 
sites have uniform strengths and are placed in the sample 
at random, but non-overlapping, positions. The pinning 
strength fp is varied from 0.2/o to l.O/o, and dp is varied 
from A/3 to A (i.e., the pin density Up varies from l/A^ 
to 9/A2). 



III. MAGNETIC FLUX DENSITY PROFILES 

Several general features of our simulations are shown 
in Fig. 1. In the upper frame of Fig. la, we show a top 
view of the vortex positions after the external field has 
been ramped up from zero. As we have stated, this ex- 
ternal field is represented by the vortices in the unpinned 
regions to the left and right of the central, pinned, sample 
region. Here, vortices have been added to the unpinned 
region to a final density of about 1.2 vortices/ A^; since 
each vortex carries a flux $0, this corresponds to a mag- 
netic field of 1.2 $o/A^. For a real superconductor with 
a penetration depth of, e.g., lOOOA, this corresponds to 
H = 2.5 kOe. 

We note in Fig. 1(a) that many of the vortices added 
to the unpinned region have been forced into the central 
sample region at this stage. They do not do so uniformly 
due to the presence of 3456 pinning sites (not shown). 



2 



with a typical intersite distance of A/2 and fp = 0.9/o. 
We see the characteristic density gradient determined 
by a balancing of the vortex-vortex forces with the lo- 
cal pinning forces. Since this gradient was achieved in 
our simulation solely by the slow ramping of an exter- 
nal magnetic field, we have obtained the field profiles 
inside a pinned superconductor using only microscopic 
information such as vortex-vortex and vortex-pin inter- 
actions. We should also contrast our simulations with 
those modeling current- driven vortices. In such simula- 
tions the driving force on each vortex is somewhat artifi- 
cially modeled by an externally-imposed "uniform" cur- 
rent. Our simulation correctly models the driving force 
as a result of local interactions. 

The lower frame of Fig. la shows the resulting flux den- 
sity profiles, found by averaging the vortex density over 
slices parallel to the sample edges. Such profiles clearly 
show the essentially constant flux density in the exter- 
nal regions, and the detailed nature of the flux gradient 
within the sample. Of course, these proflles may be ob- 
tained at any value of the external fleld. Figure lb shows 
the system after the external field has been ramped down 
from a high value to zero. The small field outside the 
sample is an artifact due to the smearing of the vortex 
fields. Now, flux remains trapped within the sample and 
the fleld gradient has changed sign. We notice that near 
the sample edges, where the fleld is small, the gradient 
in the flux density is quite large. Thus our simulation 
correctly models the increase in flux gradient (or, cquiv- 
alently, critical current) at low fields, where intervortex 
interactions are weak and pinning dominates. 

In Fig. 2 we show fiux density profiles for a complete 
cycle of the field, with the same sample parameters as 
in Fig. 1. During the initial ramp-up stage (Fig. 2, left), 
we increase the external field from zero to a final value 
of about 1.9 ^o/X^. We see the evolution of the internal 
fiux profile from first penetration at low fields, to the flrst 
complete penetration at a fleld H* « 0.8 $o/A^, to higher 
values of B at larger H. We again note the flux gradient 
is quite high at low flelds, but becomes flatter — and less 
field-dependent — at high fields. 

Of course, in real superconductors no vortices will en- 
ter the sample until H > Hd « (lnK/47r)(<I>o/A^), where 
K = X/^. However, for k's in the wide physically rele- 
vant range from 2 to 100, Hd varies from 0.05 ^o/X^ to 
0.36 <l>o/A^. Thus, Hd is small in the range of fields we 
explore. In any event, since we are only interested in the 
mixed state and not the Meissner phase, we will work in 
the approximation where Hd is negligible. 

During the ramp-down stage (Fig. 2, center), the field 
is lowered through zero to large negative values. The 
ramping down is initially effected by simply removing 
vortices from the unpinned region. However, after the 
external fleld reaches zero, it is reversed by the addition 
of antivortices in the unpinned region. During the begin- 
ning of this ramp-down stage, we note the appearance of 



the characteristic "gull- wing" flux profile as the internal 
remnant flux located close to the sample edges begins to 
be removed. Notice that at external fields near zero the 
internal fleld hardly changes at all as the external field 
is swept. This is again because of the very steep gradi- 
ents possible near zero field, where pinning dominates. 
Thus, the effect of a change in an external field near zero 
propagates only a very small distance into the sample. 

As the field decreases below H ^ (in Fig. 2, center), 
B{x) continues to have its A-shaped profile. We note that 
for small negative flelds the sample contains both vortices 
and antivortices. However, the pinning for both types is 
attractive, and so they remain locally trapped and anni- 
hilate only when their mutual attraction overcomes the 
pinning. This only occurs when they are closely spaced, 
within 0.3A. Finally, in the last ramp-up stage (Fig. 2, 
right), the full cycle is completed by increasing the field 
from the large negative value up to a large positive fleld, 
where the flux proflle looks identical to the initial ramp- 
up stage of the cycle. 

One clear advantage of our simulation is that we can 
obtain direct spatio-temporal information on the distribu- 
tion of flux inside the sample. However, experimentally 
this is quite difficult, especially for bulk samples. In- 
stead, average quantities, like magnetization curves, are 
typically obtained. From the fleld cycles shown in Fig. 2, 
we can easily obtain such magnetization loops from our 
simulation. Further, in our simulation it is simple to vary 
microscopic parameters such as pin density and strength. 
Thus, our simulations allow for a systematic study of 
the dependence of macroscopic measurements, such as 
the magnetization, on microscopic system parameters. It 
may also be possible to use our results in the reverse 
problem, so that some understanding of the microscopies 
of the pinning may be obtained from experimentally 
determined macroscopic measurements. 

IV. MAGNETIZATION HYSTERESIS LOOPS 

Experimentally, what is typically measured is the av- 
erage magnetization over the sample volume. In our sim- 
ulation, we thus calculate the average magnetization 

M=^j{H-B)dV . (5) 

In Fig. 3 we construct magnetization loops as two key 
sample microscopic parameters — the pinning density and 
strength — are varied. Fig. 3a shows complete magneti- 
zation loops obtained with the density of pins held con- 
stant at 4/A'^, but at three different values of the pin- 
ning strength fp. One can see clearly that by increasing 
the pinning strength the hysteresis loops become much 
wider. This is because a large pinning force yields a large 
fleld gradient. Thus M, which is essentially the difference 
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between the internal and external fields, will be larger 
for large fp. For instance, the remnant M is larger for 
stronger pinning. The M{H) loops all show a maximum 
when the external field is small {H < H*) and close to 
H* . This again is due to the pinning being most effective 
for low fields {H < H*). Figure 3b shows magnetization 
loops obtained for several pinning densities. Experimen- 
tally, one may systematically vary this parameter by the 
introduction of columnar defects using irradiation ||l|,|lO| . 

V. CRITICAL CURRENT VERSUS PINNING 
DENSITY AND STRENGTH 

Although magnetization loops are very useful for com- 
parison with experimental data, we have emphasized that 
our simulations allow us to directly compute the local 
flux distribution inside the sample. Thus, we may di- 
rectly measure the local critical current density Jc using 
Maxwell's equation dB/dx = fioJ. At every point on 
flux density profiles such as Fig. 2 we may compute the 
local slope (= dB/dx) and the corresponding local field 
B. This allows us to determine a large number of values 
of Jc{B). We then bin these values to obtain suitably 
averaged curves of Jc vs. B. 

As we have discussed, there are in the literature a 
great variety of functional dependences of Jc on B, cor- 
responding to different ad hoc electrodynamical assump- 
tions. The original Bean model predicts Jc to be indepen- 
dent of B. The varying slopes of the flux density in Fig. 2 
show that this prediction is not borne out in our simu- 
lation (except at relatively high-fields where the vortex- 
vortex force dominates; e.g., for weak-pinning samples 
with X^Up = 4.0, fp = 0.2/o). Kim et al. |^ have pro- 
posed that the critical current depends on B as 

a = Jc(B + 6o), (6) 

where a is field-independent and has units of force per 
unit volume. In this model, plots of 1/ Jc vs. B should 
appear as straight lines with slopes 1/a and intercept 
bo /a. The physical interpretation of the constant bo in 
Kim's model is unclear 

In Fig. 4 we plot 1/ Jc vs. B, with Jc determined from 
our flux density plots during the initial ramp up phase. 
We plot 1/ Jc for several realizations of the pinning den- 
sity Up and strength fp. Fig. 4a shows 1/Jc vs B for 
four different field sweeps with the pinning density var- 
ied from 1.0/ to 9.0/A^; in Fig. 4b we vary the pinning 
strength from 0.2/o to 0.9/o. Over a large region of the 
field, we find that 1/Jc is indeed linear in field, as in 
Kim's model. We can then fit the linear portions of each 
curve to straight lines as shown, and extract the inverse 
slope a. For fields such that B ^ bo, Kim's relation 
reads a w JcB which is the Lorentz force per unit vol- 
ume. Since this force is exactly balanced by the pinning 



force, we can interpret a as the maximum pinning force 
per unit volume, bo is typically in the range of 0.4 to 0.7 
$0/^^! but even below bo, a is clearly a measure of the 
relative effectiveness of the pinning. 

In the inset to Fig. 4a, we plot the values of a deter- 
mined from the slopes of the 1/Jc curves as a function of 
the pinning strength fp or density Up. The pinning force 
per unit volume has an approximate linear rise with Up, 
and the curve with dark triangles follows a ^ fp^ (if 
we assume that a = when fp = 0). Even though the 
vortex dynamics in our samples is not dominated by elas- 
tic flow and collective weak-pinning, it is interesting to 
compare these results with the predictions of the Larkin- 
Ovchinnikov (LO) |]ll| collective-pinning theory — where 
weakly-pinned vortices interact elastically inside a typi- 
cal correlated volume. The 2D LO prediction for rigid 
vortices becomes 

JcB ^ npfl , (7) 

which is somewhat different from 

JcB ^ npf^-^ , (8) 

obtained from our strongly pinned vortices. The oppo- 
site regime of the LO weakly-pinned collective vortices is 
given by the very strongly-pinned independent vortices 
where 

JcB ^npfl . (9) 

Thus, our results indicate that our vortices are in an 
intermediate state between the two extreme regimes de- 
scribed above. 

We plot our values for Jc in practical SI units. The 
weakest pinning in our simulation occurs at our highest 
fields, where 1/Jc is about 100^oA^/<i>o- For a A of 1000 
A, this corresponds to a critical current Jc = 1.6 x 10^ 
A/cm^, which is in practice a very reasonable value. 
Our highest critical currents, at low fields and high pin 
strength or density, are about a factor of ten higher. 
Thus, our parameters generally appear to model realistic 
materials. 

VI. CONCLUSIONS 

To summarize, we have perfomed molecular-dynamics 
simulations of vortices and antivortices interacting with 
a controlled range of pinning strengths and densities. In 
these simulations we have only considered vortex-vortex 
and vortex-pin interactions; no extra force was needed 
to simulate a Lorentz force. Thus, our results show that 
the Lorentz force can be considered as a consequence 
of a flux gradient arising strictly from the interactions 
of vortices and pins. We compute the flux density pro- 
file that develops with a varying applied field, for both 
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vortices and anti-vortices as the external field is cycled 
through a loop. Our computed complete hysteresis loops 
show realistic behaviour with varying pinning strength 
and density, indicating that our model contains the es- 
sential physics. We have obtained Jc{H) by focusing on 
the flux gradient that develops naturally from the vortex- 
pin interactions and find that it monotonically decreases 
with an increasing external field with the fall off deter- 
mined by the microscopic pinning parameters. 
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FIG. 1. Top view of the region where flux lines, indicated 
by dots, move, (a) Snapshot during the initial ramp-up 
phase, (b) snapshot of the remnant magnetization after 
ramping down the external field. The bottom panels show 
B{x) = {36X)~^ J^^^ A/ B{x,y), i.e., the flux density profile 
versus x, averaged over the vertical direction y. The 24A x 36A 
sample has 3456 pinning sites, and fp = 0.9/o. 

FIG. 2. Magnetic flux density profiles B{x) for the (1) ini- 
tial ramp-up phase, (2) ramp-down stage reaching a negative 
field, and (3) final ramp-up phase, for the same sample de- 
scribed in Fig. 1 and the text. The flat plateaus on either 
side of the sample show the density in the unpinned region, 
mimicking the external field, and the jagged V- and A-shaped 
profiles correspond to the flux density in the pinned region. 
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FIG. 3. Magnetization hysteresis curves M{H{t)). In (a) 
the maximum pinning force is varied {fp = 0.9/o, 0.55/(), 
0.2/o) for a fixed average distance between randomly dis- 
tributed pinning sites, dp = A/2 (i.e., X^Up = 4). In (b) 
the pinning-site density Up is varied while fp = 0.55/o. A 
higher value of fp and/or rip increase Jc (~ width of the 
M{H) hysteresis loop) in the manner shown in Fig. 4. For 
each M{H) loop shown, tiie maximum number of flux lines 
inside the pinned sample is about 1000. 

FIG. 4. (a) l/Jc{B) for several values of the pin density Hp 
(and fixed pinning-site strength fp); (b) l/Jc{B) for several 
values of fp (and fixed Up). The insets show the dependence 
of the maximum pinning force a on fp (dark triangles) and 
on Up (open circles) The values of a arc obtained from the 
(solid line) linear fits shown in the larger panels. 
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